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SUMMARY 



In 1970, Yefimov proposed a linearized model for predicting the effect 
of water waves on the adjacent wind flow. We investigate the impact of 
the non-linear terms neglected in his model. In this investigation, the 
sea surface is described by a single harmonic in x and the wave induced 
motion in the atmosphere is restricted to this harmonic. Buoyancy effects 
are neglected and the turbulent fluxes of momentum are approximated in terms 
of an eddy viscosity. The equations are integrated numerically until a 
steady-state is achieved. These solutions show that the nonlinear terms 
which were neglected by Yefimov are important, and that they generally 
lead to a smaller response. The numerical solutions for the mean wind are 
compared with observations obtained by Davidson and Brutsaert. Although 
the numerical solutions show some local agreement with the observations, 
the general behavior is at variance with the observations. 
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1. Introduction 



In this paper a numerical model is formulated to predict the effect of 
water waves on the adjacent air flow* Data described by Davidson (1974) and 
Brutsaert (1973) show the effect of water waves on the mean wind. Many 
studies of the air flow over waves have been carried out for the purpose of 
determining wave growth [Miles (1957), Phillips (1957), Davis (1970, 1972), 

Lee (1972)]. Reynolds and Hussain (1972) and Davis (1970, 1972), and Yefimov 
(1970) emphasized the importance of including the wave-related turbulent 
Reynolds stresses; these stresses were neglected in the linear model by 
Miles (1957). We will follow the formulation of Yefimov (1970) except that 
we will solve the initial value problem rather than the steady state equations. 
This is because we wish to include the nonlinear interaction between the wave 
induced motions and the mean flow. We shall see that these interactions are 
important . 

This study treats a single harmonic, fully developed (not affected by 
the air motion), water wave field. The wave induced motion in the atmosphere 
is also restricted to one harmonic in x , but the amplitude and phase of the 
harmonic are free to vary. Actually as shown by Lee (1972), nonlinear 
effects produce higher harmonics. We neglect these effects because we are 
mainly concerned with the more easily observed effects on the mean wind. The 
equations for a homogeneous fluid are written in terms of the Fourier coeffi- 
cients of the wave, finite differences introduced, and the equations integrated 
until a steady state is reached. These integrations are carried out for a 
number of different cases and the results are compared with the observations 
of Davidson (1974). In order to obtain better agreement with the observations, 
other boundary conditions and other diffusion coefficients are suggested. 

The nonlinear effects on the mean wind are investigated by comparison with 
the linear solutions. 
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Observational results used in comparisons were obtained from eddy correla- 
tion measurements of the momentum flux along with measurements of the mean 
wind at several levels and wave height measurements. The results were inter- 
preted (Davidson, 1974) with respect to the dependence of the drag coefficient 

on the non-dimensional parameter, c/U^ , where c corresponds to the phase 

1/2 

speed at the swell frequency and is the friction velocity, = (j/c) 

This dependence of the drag coefficient on the wave influence was extended 
in order to estimate the wave influence on the mean wind profile. The latter 
is what we examine in this study. 

2. Basic equations 

The equations for a homogeneous fluid are 



where 



dV 

V • W = - Vp + F , 

at PJ PJ 

V • V = 0 , 



F x ■ s ( - u ' 2) + & <- • 



Bx 



F = >rr (“ u'w') + (” w' 2 ) . 



( 1 ) 

( 2 ) 

(3) 

(4) 



In this formulation, the turbulent Reynolds stresses are assumed to be inde- 
pendent of y and p is the non-hydrostatic pressure divided by the 
density. The overbars represent local time averages, and u' and w' are 
the turbulent velocity fluctuations. The portion of the velocity field which 
is averaged over the turbulence time scale is given by 

V = [U(z,t) + u(x,z,t)] i + w(x,z,t) j , (5) 
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where U is the spatially averaged (mean) wind, and u and w are departures. 
Thus 



where ( ) is the average over one wavelength in x . 

We shall follow Yefimov (1970) in representing the turbulent fluxes by 



j T „■ ,r du . .chl , Bw. 

' u “ = t = k S + v( 5I + S ) ’ 



,2 ,2 

u = w 



( 6 ) 

(7) 



Note that a mathematical treatment given later by Davis (1974) tends to 
justify the eddy viscosity form (6). 

If (5), (6) and (7) are now substituted into (1), (3) and (4), and the 

curl (^ • Vx) of (1) is then taken, the following equation for the vorticity 

, du 5w. 

( § ariS6S: 



& 

at 




+ 





-n3 

o u 



,c)u , Sw., , 

+ V fe + ^ )] " V( -, 2 



Sx Sz 



a 3 w 
+ £L -|) 
a x 3 



( 8 ) 



(Note in this derivation, Su/St has been neglected, and the continuity 
equation (2) has also been used.) 

Equation (8) can further be separated into spatially stationary and non- 
stationary parts by averaging over one wavelength, then subtracting the 
averaged equation and neglecting the non-stationary perturbation products. 
This yields, for the stationary part 




+ 



4 - 



oz 



(9) 
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and, for the non- stationary * 



d>£ 

t 



+ U^I + 

OX 







3 3 

r 3 U .6 W-, 

v[ — x- - + — o] 

dx dz dx 



( 10 ) 



Since the motion is non-divergent (2), a stream function can be defined 
such that 



cH/ dilj 

u “ dz ’ W “ “ dx • 



( 11 ) 



When these relations are introduced into (10), it becomes 

♦ ■&><& + - g fi - 4n< 4 - 4)i- '-[-4* - 4 



( 12 ) 



Now we introduce sinusoidal variations and look for solutions which move 
with the speed c of the water waves 



ijr = A(z,t) cos [k(x-ct)] + B(z,t) sin [k(x-ct)] . (13) 



Substitute (13) into (12) and separate the sine and cosine terms. Equating 
the coefficients of the cosine terms gives 



- k 2 A] = k(U-c) (Bk 2 

dz 



d 2 B x , . d 2 U d 2 r ,d 2 A , 2 , 2 f d 2 A . 2 1 

— =■)+ Bk — 7 : + — ~[V ( — -z + Ak )]+ vk [ — x + Ak ] . 

dz dz dz dz dz 



(14) 



The sine terms give 



d r d 2 B , 2„, , w ., 2 d 2 A, .,d 2 U . d 2 r ..,d 2 B 2, , 2 r d 2 B . „,2, 

Vj [ — 9 “k B]= - k(U-c) (Ak r)-Ak — - + — r[V( — r- +Bk )]+Vk [ — « + Bk ] . (15) 

dz dz dz dz dz dz 
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When (11) and (13) are inserted into (9) we obtain 




( 16 ) 



Equations (14) - (16) clearly represent a coupled, non-linear system, in 
S 2 U 

that terms such as B occur. This essentially precludes analytical solu- 



tion to the full set of equations, Yefimov proceeded by first considering the 
functions involved to be independent of time, thus reducing the system to one 
of ordinary differential equations for the steady state behavior. Furthermore, 
in (14) - (15) he approximated the mean flow [our U(z,t) term] by the undisturbed log 
profile, (U 5 V /h) In (z/z q ) , thus arriving at a linear system. After solving 
this system, by what he terms a ’’matrix distillat ion’ 1 method, he then solved 
for the deviation from the log profile by direct integration of the steady-state 
counterpart of (16). As noted in the introduction, a major part of our investi- 
gation concerns the effect of these neglected non-linear terms. Also note, as 
we shall return to it later, that (16) can be integrated twice to give a simpler 

„ . - Su 

equation for ^ . 

3. Boundary Conditions 

Determination of the correct boundary conditions for (14) - (16), specifically 
at the lower boundary is not as straightforward as the derivation of the 
equations. Yefimov (1970) defined the boundary conditions on the free surface 
at the lower boundary to be 




U (T] , t ) + u(x,T],t) = - akc cos k(x-ct) , 



(17) 



w(x,r|,t) = akc sin k(x-ct) , 



where the height of the free surface is given by 



T| = a cos k(x-ct) . 



(18) 
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1 / 2 

Here a is the wave amplitude and c = (g/k) . He then assumed that 

the values of U, u and w on the free surface can be adequately approxi- 
mated by their values at z = 0 , if r| is small. With u and w replaced 
by their values in terms of (11) and (13), this reduced to 



Sa 

Sz 


(0,t) = - akc , 


A(0,t) = ac 


Sb 

9z 


(0,t) = 0 , 


B(0,t) = 0 . 



However, note that (17) can be rewritten 

U(r],t) + u(x,r],t) = - k cT) , 

hence, it seems more correct that u(x,r],t) can be approximated by its 
value at z = 0 only if U(r],t) can also be (which is what Yefimov did). 
However, the Taylor expansion for U(r),t) is 

u (t) » t ) = u(0,t) + (g) o n . (20) 

Since, for the standard logarithmic profile, (§^) Q = z 9 it is very 
possible that u(x,r],t) ^ u(x,0,t) . Note that a similar difficulty with 
w(x,0,t) is not likely, since U(t],t) does not appear in this condi- 
tion. 

Since our main intent is the investigation of the effect of the non- 
linear terms in (14)- (16), we shall use (19) as our boundary condition. 

However, since our results are at variance with the observational studies, 
further consideration of alternative formulations to (19) seems warranted. 

Note, a consequence of Yefimov's boundary conditions is that, near z = 0 , 

A(z,t) = 0(1) and B(z,t) = O(z^) . 



9 



Then, if one assumes that U(z,t) is also 0(1) , as well as second order 
and higher derivatives of A and B , then (15) implies that 




, near z = 0 . 




This apparently induced Yefimov to use —z =0 at the lower boundary in 



his numerical scheme. 

Because of the requirement that the perturbations die out at the upper 
boundary (z = H) , both A(z,t) and B(z,t) , as well as their first 
derivatives are set to zero at z = H . 

4. Representation of the Mean Wave Field 

We remark that in the case of no wave (i.e., a = 0), the solution to 
(14) - (16), plus the boundary conditions, is the log profile 



Note, this representation for the log profile does not generate a singu- 
larity at z = 0 . Also note, from (16), this immediately implies a form 
for the diffusion coefficient 



Now, in the general case (a ^ 0) , there will be a departure from 
the log profile which we represent 



A(z,t) = B(z,t) = 0 



U(z,t) = X In ( 7" + D 



O 



K = K u* (z + z o ) 



( 21 ) 



U(z,t) = -£■ In (— + 1) + u(z,t) 






( 22 ) 



o 
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where 



U Vc = friction velocity, z q = surface roughness, h= von Karmon's constant. 

When (22) is substituted in equation (16), and the resulting relation integrated, 
we have 

z 

i \ <** < 23 > 

o 

where the wave induced stress is 

rr 1 1 /^B A _* /0/ v 

T w = " UW = " 2 k A “ SI B) * (24) 

We note that Yefimov (1970) replaced the diffusion coefficient K in (23) by 
a constant V , which is apparently the eddy viscosity coefficient used 
in (6). Whether this was a misprint in the translation is not clear, and 
we investigated the solutions for both K as given by (21) and with K = v • 

5. Numerical Solution Procedure 

As noted above, Yefimov solved the linearized steady-state equivalents of 
(14) - (15) by a procedure described only by the term "matrix distillation." 

In contrast, we chose to solve the time-dependent non-linear system (14), (15), 
(22) - (24) by difference equations, and time iterate until a numerical steady 
state evolved. To implement this, (14) and (15) were first reduced to difference 
equation form, using central finite differences on all spatial operations and all 
time operations except the diffusion terms (i.e., those with coefficient v )• 

The latter were evaluated at the previous time step in order to eliminate a 
diffusive instability. The left hand sides of these equations were then solved 
for the true tendencies with the Gauss elimination technique described 
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in Richtmyer (1957) . Leapfrog time differencing was used except that the 
integration was restarted every 50 time steps with the finite difference 
scheme developed by Matsuno (1966). This procedure eliminates solution 



separation. 

The values of A(z,t) and B(z,t) thus computed are used to update the 
mean wind, as given by (21) - (24) , and this updated mean wind is then used 
to compute A(z,t) and B(z,t) at the next time step. The boundary condi- 
tions take the following finite difference form 



A(0,t) = ac , 


(25) 


B(0,t) = 0 , 


(26) 


A(Az,tJ). .-..A CP.it,). = _ akc 

Az * 


(27) 


B(Az,t) ■ B(0,t) = 
Az 


(28) 


A(H,t) = B(H,t) = 0 , 


(29) 


A(H,t)- A(H-Az,t) B(H,t)- B(H-At) n 

Az Az * 


(30) 



where Az is the grid increment. 

The following initial conditions were used to start the solution 

A(z,0) = ac sinh [k(H^- z)]/sinh kH^ , (31) 

B (z , 0) = 0 , 

u * 

U(z,0) = — In (f- + 1) , (32) 

o 

where = H-Az . Note, the fields (31) represent potential flow over the 
waves . 

From these initial conditions, as noted above, the time integration of 
the system of equations is accomplished by first predicting A and B from 
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(14) and (15) and then by obtaining U from (21) - (24). The integration is 



continued until the quantity 



_H_ 

Az 



S 



T 



[ s 

i=l 



1 

[T(z.,t) 



T(z i ,t - At)] 



V /2 , 



(33) 



satisfies the condition 

-5 2 -1 

S T < 10 m sec . 



(34) 



Condition (34) is taken to describe steady state, 

6. Observational studies 

Our numerical results will be compared to the results of several observa- 
tional studies which have considered this problem and have looked for a 
dependence on c/U^ . Brutsaert (1973) has suggested that the wind shear 
should have the following dependence 



dU 

dz 



U * c 2, 

= ^ [1+ ^“ a) 1 



(35) 



for c/U^. > a • From the data he found that 

(3 = .006 and a = 29 . 

Davidson (1974) examined the wind at 10 meters with the BOMEX data and 
obtained the following relation 

U, 

U = TJ [ In (~~) + 6.57 l + .164 - 26.3)] . (36) 

o * 

These studies show increased wind due to the wave for c/U^. > CL where a is 
between 25 and 30. We will see later that the assumption that ii is dependent 
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on only c/U^ is probably too restrictive. In particular Brutsaert (1973) 
assumed that the wave stress is constant in height which is neither reasonable 
nor consistent with observations. 

7. Numerical results 

The numerical model, as described above, was integrated for a number of 
the waves considered by Yefimov (1970) and, as a test, the system was made 
linear by setting u equal to zero. When the dimensions were restored to 
Yefimov’s solution, the comparison was quite good. It should be noted here 
that Yefimov non-dimens ionalized his variables by scaling, and the results he 
presented corresponded to unsealed values of U* in the range 0.01 msec ^ < 

< 0.05 msec \ or approximately an order of magnitude smaller than the 
values for which Davidson and Brutsaert obtained experimental data. The results 
which will be shown use the following values for the various parameters, unless 
otherwise noted: 



AZ = 0.25m 
At = 0.0125sec 
= 20m 
K = 0.35 

z = 3 x 10 ^m (37) 

o 

a = 0.1/k 

The value for was selected after numerical experiments showed larger 

values did not significantly alter the wind field at the values of c and U^. 
of interest. 

Fig. 1 shows a typical stress field prediction for the linear (Yefimov) 
model for different values of U^, at wave number k = 0.25m ^ and 
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2 -1 

V = 0.24 m sec . The values for k and U are typical for the BOMEX 

& 

data used by Davidson, and the value for v is that used by Yefimov. All 
of the curves have the same general pattern, with initially growing ampli- 
tude with increasing U*, up to some maximum amplitude, followed by decreas- 
ing amplitude for further increases in U^. Yefimov referred to this 
as a resonance-like behavior. 

Fig. 2 presents the stress field predicted for the same wave number and 
values of using the full non-linear model. Observe the qualitative 

behaviors are similar; however, two important quantitative differences emerge. 
First, the magnitude of the maximum stress amplitude for the non-linear model 
is only about 35^ of that predicted by the linear model. We found this 
pattern consistently through all ranges of wave numbers considered, with 
the maximum stress amplitude predicted by the non-linear model usually only 
one-quarter to one-half that of the linear model. Secondly, note the apparent 
"resonance 11 is noticeably less pronounced in the non-linear case. This 
pattern also appeared consistently over the range of wave numbers. From 
these observations, we have concluded that the presence of the non-linear 
terms in (14) - (16) acts strongly to damp the wave induced stresses, and 
that to approximate (14) - (16) linearly, as Yefimov did, is not valid. 

Some further calculations were carried out with ak = .05 rather than the 
value ak = 0.1 that was used by Yefimov. With this relation for a the 
comparison between the linear and nonlinear solutions was closer except 
for the longer waves near k = .05 m ^ . As a final point, note that both 
models predict the main contributions to the stress field occur in the lowest 
5-8 meters. This is especially important, since both experimental studies 
referred to above have limited or non-existent data in this region. 
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In comparing these figures to the wind shear hypothesized by Brutsaert, 
(35), it is obvious that neither the linear nor the non-linear models agree 
with his assumption of constant stress at all heights. As noted above, we 
feel this assumption was too restrictive, though until additional data is 
obtained in the region below 8 meters it cannot be fully discounted. 

To compare the results of the numerical models with Davidson 1 s, we ex- 
amined the deviation from the logarithmic profile predicted by Davidson: 



~ tt * -i /10\ 

u io ■ u io - oi ln ( r> 

o 



&3T (c - 26 ‘ 3 u *> - 



(38) 



c > 27.5 U. 



for the neutral case, with our predicted steady state value of u(10,t) , 
over a range of representative values of c and U* . Fig. 3 shows the level 
curves of u^q predicted by Davidson. Note that, even allowing for uncer- 
tainties in the determination of numerical coefficients, Davidson predicts 
essentially linear isolines. Fig. 4 shows the isolines determined by our 
non-linear model, using (21) to determine the diffusion coefficient, K . 
Although there are obviously regions, especially near the line c = 27.5 , 

where locally linear isolines occur, it is clear that our model will not 
produce this behavior globally. More seriously, although generally our 
model predicts increasing magnitude for the deviation as (c/U*) increases, 
which is consistent with Davidson's observations, the algebraic sign is not 
consistent for observations show a positive deviation growing in magnitude, 
while our model predicts a negative deviation growing in magnitude. 

The observed sign of u at 10 meters could be obtained by reversing 
the curves in Fig. 2 or by distorting them in such a way that the negative 
area would be larger than the positive area. Davidson (1974) did observe 
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negative stresses at the 3 meter level and in an earlier study with Lake 
Michigan data Davidson and Frank (1973) found positive stresses in the 
lower layers. It is not clear at present how our model could be restruc- 
tured to reverse this algebraic sign, although we suspect that the 
boundary conditions are the major influence. There are two specific areas 
of the boundary conditions that we feel need more study. The first of 
these is the adequacy of replacing values for the horizontal components 
on the free surface by their values on the mean surface, as was discussed 
in section 3. A second point is the validity of the parameterization of 
the turbulent stresses used by Yefimov (see Eqs. (6) and (7)). The be- 
havior of these stresses near the water surface is particularly crucial. 

Figs. 5, 6, and 7 show the isolines predicted by other variants on 
the model. Those in Fig. 5 are derived from the linear model. Those in 
Fig. 6 are obtained from our non-linear model, by replacing the variable 
diffusion coefficient, K in (23), by the constant V , as Yefimov T s paper 
indicated. The curves in Fig. 7 are derived using a V which varied 
with according to 

V = 1.2 sec -1 U* . (39) 

This was suggested by the development of Davis (1974). Except for the 
linear model’s unacceptably large predictions for the deviation (due to the 
’’damping 1 ' we hypothesized in the non-linear model) , all of these share the 
same qualitative features as we discussed for Fig. 4. This essentially 
leads us to believe, as mentioned above, that any changes in algebraic sign 
of the deviations must be through changes in the boundary conditions, or 
perhaps through the formulation of the wave induced turbulent diffusion. 
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8. Conclusions 



The purpose of this paper was to investigate the linear model proposed 
by Yefimov from two aspects* First, to determine the effect of the non- 
linear terms neglected by him, and secondly to determine the agreement 
between predictions by the linear and non-linear models and observational 
studies by Davidson (1974) and Brutsaert (1973), which appeared after 
Yefimov's work. 

The main conclusions of this study are that the non-linear terms neg- 
lected by Yefimov are in fact quite important when ak = 0.1 , and impart 
a strong damping to the induced fields, that neither a linear nor a non- 
linear model will produce a height-independent stress as predicted by 
Brutsaert, and that although the non-linear model produces some agreement 
with Davidson's observations, there are also important disagreements, the 
most important being in algebraic sign. However, since the model predicts 
the major contribution to the induced stress to occur in the region near 
the mean surface, and since there are significant questions about the 
proper type of the mathematical boundary conditions, resolution of the 
differences and determination of how the model should be modified to pro- 
duce results more consistent with observation demands more complete obser- 
vational data on the induced stresses near the surface. 
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